library(dplyr)
library(raster)
library(rhdf5)
library(rgdal)
library(neonAOP)
library(reshape2)
library(ggplot2)
library(MASS)
options(stringsAsFactors = FALSE)
Import D17-SOAP Field Data
## read in file
insitu.SOAP<-read.csv("/Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/insitu/veg-structure/D17_2013_SOAP_vegStr.csv")
## visualize structure
str(insitu.SOAP)
## 'data.frame': 1231 obs. of 30 variables:
## $ siteid : chr "SOAP" "SOAP" "SOAP" "SOAP" ...
## $ sitename : chr "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" ...
## $ plotid : chr "SOAP43" "SOAP43" "SOAP43" "SOAP43" ...
## $ easting : num 297065 297070 297056 297067 297057 ...
## $ northing : num 4100726 4100722 4100722 4100719 4100729 ...
## $ taxonid : chr "ABCO" "ABCO" "ABCO" "ABCO" ...
## $ scientificname : chr "Abies concolor" "Abies concolor" "Abies concolor" "Abies concolor" ...
## $ indvidual_id : int 622 632 567 606 639 566 527 638 516 555 ...
## $ pointid : chr "center" "NE" "center" "center" ...
## $ individualdistance : num 5.4 5.4 4.1 8.5 5.7 4.1 8.3 6.9 6.4 6 ...
## $ individualazimuth : num 67.8 188.7 250.8 122.4 332.6 ...
## $ dbh : num 67.6 16.1 5.9 0.9 4 5.1 1.1 3.1 1 8.9 ...
## $ dbhheight : num 130 130 130 10 130 130 130 130 10 130 ...
## $ basalcanopydiam : num 8.2 0 0 0 0 0 0 0 0 0 ...
## $ basalcanopydiam_90deg: num 7.9 0 0 0 0 0 0 0 0 0 ...
## $ maxcanopydiam : num 0 3.3 2.9 0.4 2 2.4 1.3 2 0.5 3 ...
## $ canopydiam_90deg : num 0 3.2 2.8 0.3 1.5 1.9 1.4 1.9 0.5 2.9 ...
## $ stemheight : num 34.5 13 3.8 0.6 3.4 3 1.4 3.1 0.5 5 ...
## $ stemremarks : chr "" "" "" "" ...
## $ stemstatus : chr "" "" "" "" ...
## $ canopyform : chr "" "" "" "" ...
## $ livingcanopy : int 100 100 100 100 100 100 100 100 100 100 ...
## $ inplotcanopy : int 100 100 100 100 100 100 100 100 100 100 ...
## $ materialsampleid : chr "" "" "" "" ...
## $ dbhqf : int 0 0 0 0 0 0 0 0 0 0 ...
## $ stemmapqf : int 0 1 0 0 0 0 0 0 0 0 ...
## $ plant_group : logi NA NA NA NA NA NA ...
## $ common_name : logi NA NA NA NA NA NA ...
## $ aop_plot : logi NA NA NA NA NA NA ...
## $ unique_id : logi NA NA NA NA NA NA ...
#summarize
summary(insitu.SOAP)
## siteid sitename plotid easting
## Length:1231 Length:1231 Length:1231 Min. :297034
## Class :character Class :character Class :character 1st Qu.:297064
## Mode :character Mode :character Mode :character Median :298705
## Mean :298489
## 3rd Qu.:299809
## Max. :300153
##
## northing taxonid scientificname indvidual_id
## Min. :4100083 Length:1231 Length:1231 Min. : 1.0
## 1st Qu.:4100556 Class :character Class :character 1st Qu.: 304.5
## Median :4100844 Mode :character Mode :character Median : 612.0
## Mean :4100817 Mean : 632.2
## 3rd Qu.:4101027 3rd Qu.: 920.5
## Max. :4101489 Max. :1341.0
##
## pointid individualdistance individualazimuth dbh
## Length:1231 Min. : 0.000 Min. : 0.00 Min. : 0.00
## Class :character 1st Qu.: 4.700 1st Qu.: 66.95 1st Qu.: 1.20
## Mode :character Median : 6.900 Median :160.80 Median : 3.65
## Mean : 6.903 Mean :173.57 Mean : 10.10
## 3rd Qu.: 9.100 3rd Qu.:286.30 3rd Qu.: 13.68
## Max. :64.000 Max. :359.80 Max. :159.80
## NA's :249
## dbhheight basalcanopydiam basalcanopydiam_90deg maxcanopydiam
## Min. : 0.0 Min. : 0.000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 10.0 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.700
## Median :130.0 Median : 0.000 Median : 0.000 Median : 1.600
## Mean : 86.3 Mean : 7.556 Mean : 5.322 Mean : 2.355
## 3rd Qu.:130.0 3rd Qu.: 4.000 3rd Qu.: 2.000 3rd Qu.: 3.100
## Max. :140.0 Max. :220.000 Max. :93.000 Max. :80.000
## NA's :1 NA's :1
## canopydiam_90deg stemheight stemremarks stemstatus
## Min. : 0.000 Min. : -1.600 Length:1231 Length:1231
## 1st Qu.: 0.500 1st Qu.: 1.000 Class :character Class :character
## Median : 1.300 Median : 1.800 Mode :character Mode :character
## Mean : 1.863 Mean : 4.786
## 3rd Qu.: 2.400 3rd Qu.: 4.450
## Max. :70.000 Max. :134.000
## NA's :1
## canopyform livingcanopy inplotcanopy materialsampleid
## Length:1231 Min. : 0.00 Min. : 40.00 Length:1231
## Class :character 1st Qu.: 60.00 1st Qu.:100.00 Class :character
## Mode :character Median :100.00 Median :100.00 Mode :character
## Mean : 74.25 Mean : 99.78
## 3rd Qu.:100.00 3rd Qu.:100.00
## Max. :100.00 Max. :100.00
## NA's :85 NA's :84
## dbhqf stemmapqf plant_group common_name
## Min. :0 Min. :0.00000 Mode:logical Mode:logical
## 1st Qu.:0 1st Qu.:0.00000 NA's:1231 NA's:1231
## Median :0 Median :0.00000
## Mean :0 Mean :0.09748
## 3rd Qu.:0 3rd Qu.:0.00000
## Max. :0 Max. :1.00000
## NA's :1
## aop_plot unique_id
## Mode:logical Mode:logical
## NA's:1231 NA's:1231
##
##
##
##
##
Exploratory Analysis: Subset variables of interest and identify data type (Cat & Cont)
# selected potential predictors of interest at SOAP 2013
# plotid, taxonid, scientificname, dbh, dbhheight, basalCanopyDiam, Maximum canopy diameter, stemheight
# canopyform, livingcanopy, inplotcanopy
# subset dataframe-- make it not dirty
insitu.SOAP<-subset(insitu.SOAP, select=c(plotid, easting, northing, taxonid, scientificname, dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, canopyform, livingcanopy, inplotcanopy))
## subset categorical variables
SOAP.cat<-subset(insitu.SOAP, select=c(plotid, taxonid, scientificname, canopyform))
## subset continuous variables
SOAP.cont<- subset(insitu.SOAP, select=c(dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy))
Exploratory Analysis for SOAP- Frequency Categorical Statistics
## frequency summaries by:
# plotid
SOAP.cat %>%
group_by(plotid) %>%
summarise (n = n()) %>%
mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [18 x 3]
##
## plotid n rel.freq
## (chr) (int) (chr)
## 1 SOAP1343 24 1.9%
## 2 SOAP139 122 9.9%
## 3 SOAP143 114 9.3%
## 4 SOAP1515 26 2.1%
## 5 SOAP1563 14 1.1%
## 6 SOAP1611 10 0.8%
## 7 SOAP1695 106 8.6%
## 8 SOAP187 95 7.7%
## 9 SOAP223 43 3.5%
## 10 SOAP255 104 8.4%
## 11 SOAP283 14 1.1%
## 12 SOAP299 26 2.1%
## 13 SOAP331 132 10.7%
## 14 SOAP43 100 8.1%
## 15 SOAP555 6 0.5%
## 16 SOAP63 122 9.9%
## 17 SOAP95 106 8.6%
## 18 SOAP991 67 5.4%
# dominant(?) plant type
SOAP.cat %>%
group_by(scientificname) %>%
summarise (n = n()) %>%
mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [16 x 3]
##
## scientificname n rel.freq
## (chr) (int) (chr)
## 1 Abies concolor 11 0.9%
## 2 Amelanchier utahensis 2 0.2%
## 3 Arctostaphylos viscida 179 14.5%
## 4 Calocedrus decurrens 435 35.3%
## 5 Ceanothus cuneatus 49 4%
## 6 Ceanothus integerrimus 42 3.4%
## 7 Cercocarpus montanus var. glaber 40 3.2%
## 8 Corylus cornuta 3 0.2%
## 9 Pinus lambertiana 16 1.3%
## 10 Pinus ponderosa 213 17.3%
## 11 Prunus emarginata 1 0.1%
## 12 Quercus chrysolepis 141 11.5%
## 13 Quercus garryana 3 0.2%
## 14 Quercus kelloggii 65 5.3%
## 15 Rhamnus ilicifolia 12 1%
## 16 Ribes roezlii 19 1.5%
# canopyform
SOAP.cat %>%
group_by(canopyform) %>%
summarise (n = n()) %>%
mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [6 x 3]
##
## canopyform n rel.freq
## (chr) (int) (chr)
## 1 844 68.6%
## 2 Cone 155 12.6%
## 3 Cylinder 14 1.1%
## 4 Hemisphere 70 5.7%
## 5 Sphere 64 5.2%
## 6 NA 84 6.8%
Exploratory Analysis for SOAP- Summary Continuous Statistics
## Continuous variable summaries- Min/Max, Quants
summary(SOAP.cont)
## dbh dbhheight basalcanopydiam maxcanopydiam
## Min. : 0.00 Min. : 0.0 Min. : 0.000 Min. : 0.000
## 1st Qu.: 1.20 1st Qu.: 10.0 1st Qu.: 0.000 1st Qu.: 0.700
## Median : 3.65 Median :130.0 Median : 0.000 Median : 1.600
## Mean : 10.10 Mean : 86.3 Mean : 7.556 Mean : 2.355
## 3rd Qu.: 13.68 3rd Qu.:130.0 3rd Qu.: 4.000 3rd Qu.: 3.100
## Max. :159.80 Max. :140.0 Max. :220.000 Max. :80.000
## NA's :249 NA's :1 NA's :1
## stemheight livingcanopy inplotcanopy
## Min. : -1.600 Min. : 0.00 Min. : 40.00
## 1st Qu.: 1.000 1st Qu.: 60.00 1st Qu.:100.00
## Median : 1.800 Median :100.00 Median :100.00
## Mean : 4.786 Mean : 74.25 Mean : 99.78
## 3rd Qu.: 4.450 3rd Qu.:100.00 3rd Qu.:100.00
## Max. :134.000 Max. :100.00 Max. :100.00
## NA's :85 NA's :84
## Continuous variable summaries- Std Devs
sapply(SOAP.cont, sd, na.rm=TRUE)
## dbh dbhheight basalcanopydiam maxcanopydiam
## 14.927547 58.329815 19.304295 3.803459
## stemheight livingcanopy inplotcanopy
## 8.526687 41.145058 2.863307
# matrix scatterplot
plot(SOAP.cont)
#side by side boxplots
boxplot(SOAP.cont, las=2)
Exploratory Analysis for SOAP by Sampling Plots
## by species types by PlotID
xtabs(~plotid+scientificname, data=SOAP.cat)
## scientificname
## plotid Abies concolor Amelanchier utahensis Arctostaphylos viscida
## SOAP1343 0 0 12
## SOAP139 0 0 14
## SOAP143 0 0 24
## SOAP1515 0 0 23
## SOAP1563 0 0 14
## SOAP1611 0 0 10
## SOAP1695 0 0 36
## SOAP187 0 0 0
## SOAP223 0 0 8
## SOAP255 0 0 15
## SOAP283 0 0 7
## SOAP299 0 0 3
## SOAP331 4 2 0
## SOAP43 7 0 0
## SOAP555 0 0 6
## SOAP63 0 0 6
## SOAP95 0 0 0
## SOAP991 0 0 1
## scientificname
## plotid Calocedrus decurrens Ceanothus cuneatus Ceanothus integerrimus
## SOAP1343 1 5 0
## SOAP139 72 0 1
## SOAP143 13 0 34
## SOAP1515 0 0 1
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 0 1 0
## SOAP187 81 0 0
## SOAP223 0 2 0
## SOAP255 14 0 0
## SOAP283 3 0 1
## SOAP299 3 0 0
## SOAP331 85 0 0
## SOAP43 63 0 0
## SOAP555 0 0 0
## SOAP63 38 0 0
## SOAP95 61 0 5
## SOAP991 1 41 0
## scientificname
## plotid Cercocarpus montanus var. glaber Corylus cornuta
## SOAP1343 2 0
## SOAP139 0 0
## SOAP143 0 0
## SOAP1515 0 0
## SOAP1563 0 0
## SOAP1611 0 0
## SOAP1695 0 0
## SOAP187 1 0
## SOAP223 16 0
## SOAP255 0 0
## SOAP283 0 0
## SOAP299 0 0
## SOAP331 0 0
## SOAP43 0 3
## SOAP555 0 0
## SOAP63 0 0
## SOAP95 0 0
## SOAP991 21 0
## scientificname
## plotid Pinus lambertiana Pinus ponderosa Prunus emarginata
## SOAP1343 0 1 0
## SOAP139 3 6 0
## SOAP143 0 38 0
## SOAP1515 0 1 0
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 0 4 0
## SOAP187 0 8 0
## SOAP223 0 0 0
## SOAP255 0 51 0
## SOAP283 0 0 0
## SOAP299 0 16 0
## SOAP331 1 11 0
## SOAP43 9 6 0
## SOAP555 0 0 0
## SOAP63 0 53 0
## SOAP95 3 18 0
## SOAP991 0 0 1
## scientificname
## plotid Quercus chrysolepis Quercus garryana Quercus kelloggii
## SOAP1343 0 1 1
## SOAP139 19 0 2
## SOAP143 0 0 1
## SOAP1515 1 0 0
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 50 0 0
## SOAP187 2 0 2
## SOAP223 14 1 1
## SOAP255 15 0 9
## SOAP283 3 0 0
## SOAP299 4 0 0
## SOAP331 6 0 23
## SOAP43 0 0 9
## SOAP555 0 0 0
## SOAP63 10 0 15
## SOAP95 17 0 2
## SOAP991 0 1 0
## scientificname
## plotid Rhamnus ilicifolia Ribes roezlii
## SOAP1343 0 1
## SOAP139 0 5
## SOAP143 0 4
## SOAP1515 0 0
## SOAP1563 0 0
## SOAP1611 0 0
## SOAP1695 12 3
## SOAP187 0 1
## SOAP223 0 1
## SOAP255 0 0
## SOAP283 0 0
## SOAP299 0 0
## SOAP331 0 0
## SOAP43 0 3
## SOAP555 0 0
## SOAP63 0 0
## SOAP95 0 0
## SOAP991 0 1
## aggregate summary stats
# Re-instate PlotID
SOAP.cont$plotid<-insitu.SOAP$plotid
# group by
melted <- melt(SOAP.cont, id.vars=c("plotid"))
SOAP.aggregate<-summarise(group_by(melted, plotid, variable),
mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE),
min=min(value, na.rm=TRUE), max=max(value, na.rm=TRUE))
SOAP.aggregate
## Source: local data frame [126 x 6]
## Groups: plotid [?]
##
## plotid variable mean sd min max
## (chr) (fctr) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 dbh 8.400000 6.437391 1.0 15.4
## 2 SOAP1343 dbhheight 120.000000 33.879582 10.0 130.0
## 3 SOAP1343 basalcanopydiam 29.287500 30.840844 0.0 140.0
## 4 SOAP1343 maxcanopydiam 3.920833 5.736684 0.6 30.0
## 5 SOAP1343 stemheight 2.112500 1.212368 0.7 6.5
## 6 SOAP1343 livingcanopy 71.625000 33.217089 0.0 100.0
## 7 SOAP1343 inplotcanopy 100.000000 0.000000 100.0 100.0
## 8 SOAP139 dbh 7.274510 19.207526 0.4 159.8
## 9 SOAP139 dbhheight 80.983607 59.089648 10.0 130.0
## 10 SOAP139 basalcanopydiam 5.311475 13.882674 0.0 60.0
## .. ... ... ... ... ... ...
## WHOA-- That's a big spicy meatball!
## Break it down by plot again-- this time stratify by parameter
#overside columnn width restriction
options(dplyr.width = Inf)
# dbh
insitu_dbh <- SOAP.cont %>%
group_by(plotid) %>%
summarise(dbh.min = min(dbh, na.rm=TRUE), dbh.max = max(dbh, na.rm=TRUE), dbh.mean=mean(dbh, na.rm=TRUE), dbh.median=median(dbh, na.rm=TRUE),dbh.sd=sd(dbh, na.rm=TRUE))
insitu_dbh
## Source: local data frame [18 x 6]
##
## plotid dbh.min dbh.max dbh.mean dbh.median dbh.sd
## (chr) (dbl) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 1.0 15.4 8.400000 8.60 6.437391
## 2 SOAP139 0.4 159.8 7.274510 1.70 19.207526
## 3 SOAP143 0.4 44.0 6.018889 1.40 11.503276
## 4 SOAP1515 1.0 24.0 11.646154 10.50 6.357753
## 5 SOAP1563 0.0 0.0 0.000000 0.00 0.000000
## 6 SOAP1611 0.0 0.0 0.000000 0.00 0.000000
## 7 SOAP1695 0.0 60.7 7.654667 1.90 12.417656
## 8 SOAP187 1.2 102.9 8.220652 2.75 16.507677
## 9 SOAP223 0.6 107.9 16.823529 10.20 26.028387
## 10 SOAP255 0.4 78.3 10.751724 7.30 13.063271
## 11 SOAP283 0.0 33.5 6.492857 0.00 12.273578
## 12 SOAP299 1.1 64.1 32.685000 33.90 21.016767
## 13 SOAP331 0.4 72.1 10.714729 5.90 12.610971
## 14 SOAP43 0.7 84.0 11.066316 5.80 16.032312
## 15 SOAP555 0.0 0.0 0.000000 0.00 0.000000
## 16 SOAP63 0.2 50.8 14.746364 13.25 12.035816
## 17 SOAP95 0.3 61.3 10.762626 5.50 12.969474
## 18 SOAP991 1.4 10.5 6.420000 8.20 4.490212
## functionalization might be best here...but for now
## JV scripting
# dbhheight
insitu_dbhheight <- SOAP.cont %>%
group_by(plotid) %>%
summarise(dbhheight.min = min(dbhheight, na.rm=TRUE), dbhheight.max = max(dbhheight, na.rm=TRUE), dbhheight.mean=mean(dbhheight, na.rm=TRUE), dbhheight.median=median(dbhheight, na.rm=TRUE),dbhheight.sd=sd(dbhheight, na.rm=TRUE))
insitu_dbhheight
## Source: local data frame [18 x 6]
##
## plotid dbhheight.min dbhheight.max dbhheight.mean dbhheight.median
## (chr) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 10.0 130 120.00000 130.0
## 2 SOAP139 10.0 130 80.98361 130.0
## 3 SOAP143 10.0 130 59.47368 10.0
## 4 SOAP1515 1.0 130 70.88462 130.0
## 5 SOAP1563 0.0 0 0.00000 0.0
## 6 SOAP1611 0.0 0 0.00000 0.0
## 7 SOAP1695 0.0 130 66.36604 38.5
## 8 SOAP187 2.0 130 69.07368 10.0
## 9 SOAP223 10.0 130 94.18605 130.0
## 10 SOAP255 10.0 140 109.32692 130.0
## 11 SOAP283 0.0 130 28.18571 0.0
## 12 SOAP299 10.0 130 125.38462 130.0
## 13 SOAP331 10.0 130 83.65909 130.0
## 14 SOAP43 10.0 130 103.60000 130.0
## 15 SOAP555 0.0 0 0.00000 0.0
## 16 SOAP63 0.8 130 111.23607 130.0
## 17 SOAP95 10.0 130 85.84906 130.0
## 18 SOAP991 10.0 130 121.04478 130.0
## dbhheight.sd
## (dbl)
## 1 33.87958
## 2 59.08965
## 3 59.33022
## 4 65.13944
## 5 0.00000
## 6 0.00000
## 7 61.59260
## 8 60.22498
## 9 55.12896
## 10 45.67221
## 11 55.19133
## 12 23.53394
## 13 58.62428
## 14 49.95998
## 15 0.00000
## 16 43.87607
## 17 58.14379
## 18 31.77260
# basalcanopydiam
insitu_basalcanopydiam <- SOAP.cont %>%
group_by(plotid) %>%
summarise(basalcanopydiam.min = min(basalcanopydiam, na.rm=TRUE), basalcanopydiam.max = max(basalcanopydiam, na.rm=TRUE), basalcanopydiam.mean=mean(basalcanopydiam, na.rm=TRUE), basalcanopydiam.median=median(basalcanopydiam, na.rm=TRUE),basalcanopydiam.sd=sd(basalcanopydiam, na.rm=TRUE))
insitu_basalcanopydiam
## Source: local data frame [18 x 6]
##
## plotid basalcanopydiam.min basalcanopydiam.max basalcanopydiam.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.0 140 29.2875000
## 2 SOAP139 0.0 60 5.3114754
## 3 SOAP143 0.0 40 2.3771930
## 4 SOAP1515 0.0 117 23.8461538
## 5 SOAP1563 2.0 90 27.2142857
## 6 SOAP1611 1.1 220 42.6100000
## 7 SOAP1695 0.0 94 10.3000000
## 8 SOAP187 0.0 48 1.0315789
## 9 SOAP223 0.0 58 11.3255814
## 10 SOAP255 0.0 80 5.5000000
## 11 SOAP283 0.0 31 14.0000000
## 12 SOAP299 0.0 125 15.9230769
## 13 SOAP331 0.0 60 0.5174242
## 14 SOAP43 0.0 38 0.8320000
## 15 SOAP555 7.0 130 90.8333333
## 16 SOAP63 0.0 61 2.3360656
## 17 SOAP95 0.0 30 0.7809524
## 18 SOAP991 0.0 163 34.6567164
## basalcanopydiam.median basalcanopydiam.sd
## (dbl) (dbl)
## 1 22.5 30.840844
## 2 0.0 13.882674
## 3 0.0 6.012362
## 4 20.0 21.502916
## 5 24.0 22.874839
## 6 18.0 64.730286
## 7 0.0 17.087255
## 8 0.0 6.104502
## 9 6.0 15.085030
## 10 0.0 15.270379
## 11 13.0 11.948608
## 12 0.0 35.305436
## 13 0.0 5.253389
## 14 0.0 4.755664
## 15 110.0 48.803347
## 16 0.0 8.759987
## 17 0.0 3.855518
## 18 24.0 31.471754
# maxcanopydiam
insitu_maxcanopydiam <- SOAP.cont %>%
group_by(plotid) %>%
summarise(maxcanopydiam.min = min(maxcanopydiam, na.rm=TRUE), maxcanopydiam.max = max(maxcanopydiam, na.rm=TRUE), maxcanopydiam.mean=mean(maxcanopydiam, na.rm=TRUE), maxcanopydiam.median=median(maxcanopydiam, na.rm=TRUE),maxcanopydiam.sd=sd(maxcanopydiam, na.rm=TRUE))
insitu_maxcanopydiam
## Source: local data frame [18 x 6]
##
## plotid maxcanopydiam.min maxcanopydiam.max maxcanopydiam.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.6 30.0 3.920833
## 2 SOAP139 0.2 30.0 1.752459
## 3 SOAP143 0.2 11.5 1.240351
## 4 SOAP1515 0.7 5.4 2.334615
## 5 SOAP1563 0.5 3.9 1.885714
## 6 SOAP1611 1.0 5.1 2.240000
## 7 SOAP1695 0.1 60.0 3.513208
## 8 SOAP187 0.0 80.0 2.309474
## 9 SOAP223 0.3 10.1 2.576744
## 10 SOAP255 0.0 12.6 3.031731
## 11 SOAP283 0.0 3.8 1.807143
## 12 SOAP299 0.0 10.1 4.330769
## 13 SOAP331 0.0 8.5 1.934351
## 14 SOAP43 0.0 18.1 2.440000
## 15 SOAP555 1.3 4.0 2.200000
## 16 SOAP63 0.0 6.8 2.604918
## 17 SOAP95 0.0 12.4 2.096226
## 18 SOAP991 0.4 4.7 1.962687
## maxcanopydiam.median maxcanopydiam.sd
## (dbl) (dbl)
## 1 2.85 5.7366843
## 2 0.90 3.2085684
## 3 0.60 1.8008709
## 4 2.10 1.0403624
## 5 1.80 1.0021954
## 6 1.80 1.3150412
## 7 2.20 7.1522634
## 8 0.80 8.3154186
## 9 1.90 2.1445209
## 10 2.40 2.3197449
## 11 1.70 1.2256642
## 12 4.25 2.9714332
## 13 1.40 1.6596928
## 14 1.95 2.4516743
## 15 1.95 0.9423375
## 16 2.45 1.6878688
## 17 0.90 2.2582304
## 18 2.00 0.9766727
# stemheight
insitu_stemheight <- SOAP.cont %>%
group_by(plotid) %>%
summarise(stemheight.min = min(stemheight, na.rm=TRUE), stemheight.max = max(stemheight, na.rm=TRUE), stemheight.mean=mean(stemheight, na.rm=TRUE), stemheight.median=median(stemheight, na.rm=TRUE),stemheight.sd=sd(stemheight, na.rm=TRUE))
insitu_stemheight
## Source: local data frame [18 x 6]
##
## plotid stemheight.min stemheight.max stemheight.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.7 6.5 2.1125000
## 2 SOAP139 0.5 120.0 4.7459016
## 3 SOAP143 0.5 19.7 2.5307018
## 4 SOAP1515 0.6 3.4 1.7769231
## 5 SOAP1563 0.0 1.9 1.0714286
## 6 SOAP1611 0.6 4.2 1.7200000
## 7 SOAP1695 0.6 17.0 3.0566038
## 8 SOAP187 0.5 134.0 4.8789474
## 9 SOAP223 0.5 8.2 2.4441860
## 10 SOAP255 0.0 44.7 5.1019231
## 11 SOAP283 0.0 13.8 3.3285714
## 12 SOAP299 1.0 33.1 13.8961538
## 13 SOAP331 0.5 39.2 5.0098485
## 14 SOAP43 0.5 51.1 5.0470000
## 15 SOAP555 0.8 1.4 0.9666667
## 16 SOAP63 0.6 33.0 9.2327869
## 17 SOAP95 0.5 28.1 6.0509434
## 18 SOAP991 -1.6 3.8 1.8537313
## stemheight.median stemheight.sd
## (dbl) (dbl)
## 1 2.05 1.2123683
## 2 1.50 13.0803565
## 3 1.00 4.1863057
## 4 1.70 0.5867249
## 5 1.20 0.5426836
## 6 1.30 1.2318008
## 7 2.20 2.8116770
## 8 1.40 14.9537862
## 9 1.60 1.9816182
## 10 2.35 6.7463328
## 11 1.80 4.1795157
## 12 10.30 12.8098550
## 13 2.40 6.4031226
## 14 2.95 7.8090010
## 15 0.90 0.2250926
## 16 5.00 9.0170403
## 17 1.70 8.2020583
## 18 1.90 0.9241422
# livingcanopy
insitu_livingcanopy <- SOAP.cont %>%
group_by(plotid) %>%
summarise(livingcanopy.min = min(livingcanopy, na.rm=TRUE), livingcanopy.max = max(livingcanopy, na.rm=TRUE), livingcanopy.mean=mean(livingcanopy, na.rm=TRUE), livingcanopy.median=median(livingcanopy, na.rm=TRUE),livingcanopy.sd=sd(livingcanopy, na.rm=TRUE))
insitu_livingcanopy
## Source: local data frame [18 x 6]
##
## plotid livingcanopy.min livingcanopy.max livingcanopy.mean
## (chr) (int) (int) (dbl)
## 1 SOAP1343 0 100 71.62500
## 2 SOAP139 0 100 79.71311
## 3 SOAP143 0 100 95.92105
## 4 SOAP1515 20 100 82.42308
## 5 SOAP1563 NA NA NaN
## 6 SOAP1611 NA NA NaN
## 7 SOAP1695 0 100 88.57576
## 8 SOAP187 0 100 23.15789
## 9 SOAP223 0 100 79.48837
## 10 SOAP255 0 100 66.89320
## 11 SOAP283 NA NA NaN
## 12 SOAP299 0 100 84.61538
## 13 SOAP331 0 100 78.78788
## 14 SOAP43 0 100 81.90000
## 15 SOAP555 NA NA NaN
## 16 SOAP63 0 100 62.58197
## 17 SOAP95 0 100 84.90566
## 18 SOAP991 0 100 71.41791
## livingcanopy.median livingcanopy.sd
## (dbl) (dbl)
## 1 80 33.21709
## 2 100 39.42405
## 3 100 18.47703
## 4 90 20.29418
## 5 NA NaN
## 6 NA NaN
## 7 100 19.37172
## 8 0 42.40793
## 9 90 24.98225
## 10 100 46.67953
## 11 NA NaN
## 12 100 36.79465
## 13 100 41.03676
## 14 100 38.57814
## 15 NA NaN
## 16 100 48.22568
## 17 100 35.96944
## 18 80 26.73957
# inplotcanopy
insitu_inplotcanopy <- SOAP.cont %>%
group_by(plotid) %>%
summarise(inplotcanopy.min = min(inplotcanopy, na.rm=TRUE), inplotcanopy.max = max(inplotcanopy, na.rm=TRUE), inplotcanopy.mean=mean(inplotcanopy, na.rm=TRUE), inplotcanopy.median=median(inplotcanopy, na.rm=TRUE),inplotcanopy.sd=sd(inplotcanopy, na.rm=TRUE))
insitu_inplotcanopy
## Source: local data frame [18 x 6]
##
## plotid inplotcanopy.min inplotcanopy.max inplotcanopy.mean
## (chr) (int) (int) (dbl)
## 1 SOAP1343 100 100 100.00000
## 2 SOAP139 90 100 99.91803
## 3 SOAP143 100 100 100.00000
## 4 SOAP1515 100 100 100.00000
## 5 SOAP1563 NA NA NaN
## 6 SOAP1611 NA NA NaN
## 7 SOAP1695 65 100 98.03030
## 8 SOAP187 100 100 100.00000
## 9 SOAP223 100 100 100.00000
## 10 SOAP255 100 100 100.00000
## 11 SOAP283 NA NA NaN
## 12 SOAP299 100 100 100.00000
## 13 SOAP331 100 100 100.00000
## 14 SOAP43 40 100 98.90000
## 15 SOAP555 NA NA NaN
## 16 SOAP63 100 100 100.00000
## 17 SOAP95 100 100 100.00000
## 18 SOAP991 100 100 100.00000
## inplotcanopy.median inplotcanopy.sd
## (dbl) (dbl)
## 1 100 0.0000000
## 2 100 0.9053575
## 3 100 0.0000000
## 4 100 0.0000000
## 5 NA NaN
## 6 NA NaN
## 7 100 6.7867965
## 8 100 0.0000000
## 9 100 0.0000000
## 10 100 0.0000000
## 11 NA NaN
## 12 100 0.0000000
## 13 100 0.0000000
## 14 100 7.7713538
## 15 NA NaN
## 16 100 0.0000000
## 17 100 0.0000000
## 18 100 0.0000000
create a qualitative landscape assessment index (completely experimental- BE NICE!)
#subset SOAP.grand data set
#s<-subset(SOAP.grand, select = c(plot.id, easting, northing, max.dbh, max.dbhheight, max.basalcanopydiam, max.maxcanopydiam, max.stemheight))
## categorize by quartiles, use quantile values from summary and assign score
# dbh
summary(SOAP.grand$max.dbh)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 41.38 62.70 65.91 88.72 159.80
dbh.cat <- cut(SOAP.grand$max.dbh, breaks = c(-1, 41.38, 62.70, 88.72, 159.8))
dbh.cat<-unclass(dbh.cat)
#basalcanopydiam
summary(SOAP.grand$max.basalcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30.00 39.50 59.00 77.83 77.00 220.00
max.basalcanopydiam.cat<- cut(SOAP.grand$max.basalcanopydiam, breaks= c(0, 39.50, 59.0, 77.00, 222.00))
max.basalcanopydiam.cat<- unclass(max.basalcanopydiam.cat)
# maxcanopydiam
summary(SOAP.grand$max.maxcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3.80 5.95 10.10 17.51 15.25 80.00 1
## impute missing value at Plot SOAP331
SOAP.grand$max.maxcanopydiam[8] <- mean(SOAP.grand$max.maxcanopydiam, na.rm=TRUE)
# determine class quantiles scores
summary(SOAP.grand$max.maxcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.800 6.375 10.800 17.510 17.660 80.000
max.maxcanopydiam.cat<-cut(SOAP.grand$max.maxcanopydiam, breaks= c(0, 6.375, 10.80, 17.660, 80.00))
max.maxcanopydiam.cat<-unclass(max.maxcanopydiam.cat)
#stem height
summary(SOAP.grand$max.stemheight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.80 12.40 30.55 40.68 42.18 134.00
max.stemheight.cat<-cut(SOAP.grand$max.stemheight, breaks= c(0,12.40,30.55, 42.18, 134))
max.stemheight.cat<-unclass(max.stemheight.cat)
s<-data.frame(dbh.cat, max.basalcanopydiam.cat,max.basalcanopydiam.cat,max.stemheight.cat)
tree.score<-rowSums(s)
s$tree.score<-tree.score
s$tree.score.perc<- (tree.score/16)
## add on to SOAP.grand dataset
SOAP.grand<-c(SOAP.grand,s)
boxplot(s$tree.score.perc)
## simple linear regression using tree score ()
NDVI.score<-summary(lm(tree.score.perc~SOAP_NDVI, data=SOAP.grand))
## also- makesure you append CHM values from other data stack to model
## import chm
SOAP.chm <- raster("../NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif")
SOAP.chm
## class : RasterLayer
## dimensions : 1516, 3292, 4990672 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 296906, 300198, 4100038, 4101554 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif
## names : SOAP_lidarCHM
## values : 0, 68.83 (min, max)
# Insitu sampling took place within 40m x 40m square plots, so we use a 20m radius.
# Note that below will return a dataframe containing the max height
# calculated from all pixels in the buffer for each plot (NEON Website)
## we will use these values to populate in-situ data set
chm.max <- extract(SOAP.chm,
SOAP_plots,
buffer = 20,
fun=(max),
sp=TRUE,
stringsAsFactors=FALSE)
# CHM distribution visualization for each plot
# created histograms for all 14 plots of data
cent_ovrList <- extract(SOAP.chm,SOAP_plots,buffer = 20)
for (i in 1:14) {
hist(cent_ovrList[[i]], main=(paste("plot",i)))
}
## merge data into data frame, can use approach to create spatial object
## But purposes of capstone will conduct "boring regression analysis"
## merge grouped data (insitu.max) and extracted NDVI (NDVI.max)
## using plotid to merge
## first must append NDVI.max file with SOAP
#create data frame for NDVI.max
chm.max.df<-chm.max@data
# add SOAP prefix to IDs
chm.max.df$plot.id<- "SOAP"
chm.max.df$plot.id<-(sapply(chm.max.df$plot.id, paste0, chm.max.df$ID))[,1]
# rename column
colnames(chm.max.df)[colnames(chm.max.df) == 'SOAP'] <- 'Plot.ID'
SOAP.grand<-merge(chm.max.df, SOAP.grand,
by.x = 'plot.id',
by.y = 'plot.id')
# F* Yes
## Response variables: dbh, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy
## Predictor variable: SOAP_lidarCHM
## Simple Linear Fit
mod.chm<-lm(max.dbh~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.dbh, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" DBH (cm)")
abline(mod.chm)
summary(mod.chm)
##
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.517 -16.469 -8.646 0.459 78.467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9555 30.4248 0.163 0.874
## SOAP_lidarCHM 1.8172 0.8449 2.151 0.057 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.35 on 10 degrees of freedom
## Multiple R-squared: 0.3163, Adjusted R-squared: 0.2479
## F-statistic: 4.626 on 1 and 10 DF, p-value: 0.05698
mod2.chm<-lm(max.basalcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.basalcanopydiam, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Base Canopy Diameter (cm)")
abline(mod2.chm)
summary(mod2.chm)
##
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.543 -47.516 -5.578 19.236 110.410
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 133.947 46.005 2.912 0.0155 *
## SOAP_lidarCHM -1.673 1.278 -1.309 0.2197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.99 on 10 degrees of freedom
## Multiple R-squared: 0.1464, Adjusted R-squared: 0.06102
## F-statistic: 1.715 on 1 and 10 DF, p-value: 0.2197
mod3.chm<-lm(max.maxcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.maxcanopydiam, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Maximum canopy diameter(m)")
abline(mod3.chm)
summary(mod3.chm)
##
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.240 -9.587 -3.554 0.542 49.986
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.5842 14.4213 -0.803 0.4405
## SOAP_lidarCHM 0.8674 0.4005 2.166 0.0556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.18 on 10 degrees of freedom
## Multiple R-squared: 0.3193, Adjusted R-squared: 0.2512
## F-statistic: 4.691 on 1 and 10 DF, p-value: 0.05556
mod4.chm<-lm(max.stemheight~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.stemheight, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Stem Height (m)")
abline(mod4.chm)
summary(mod4.chm)
##
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.832 -17.010 -5.808 9.709 59.479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -38.0415 23.7283 -1.603 0.13997
## SOAP_lidarCHM 2.3470 0.6589 3.562 0.00516 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.91 on 10 degrees of freedom
## Multiple R-squared: 0.5592, Adjusted R-squared: 0.5151
## F-statistic: 12.69 on 1 and 10 DF, p-value: 0.005165
## add on to SOAP.grand dataset
SOAP.grand<-c(SOAP.grand,s)
boxplot(s$tree.score.perc)
## simple linear regression using tree score
chm.score<-summary(lm(tree.score.perc~SOAP_lidarCHM, data=SOAP.grand))
chm.score
##
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25200 -0.11235 0.04449 0.08071 0.24605
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.420212 0.129625 3.242 0.00884 **
## SOAP_lidarCHM 0.006105 0.003600 1.696 0.12072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1634 on 10 degrees of freedom
## Multiple R-squared: 0.2234, Adjusted R-squared: 0.1458
## F-statistic: 2.877 on 1 and 10 DF, p-value: 0.1207
## Full model in respect to dbh
mod.full.dbh <- summary(lm(max.dbh~SOAP_lidarCHM + SOAP_NDVI, data=SOAP.grand))
(mod.full.dbh)
##
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.540 -16.487 -11.476 -2.299 74.451
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -137.057 241.573 -0.567 0.584
## SOAP_lidarCHM 1.505 1.020 1.475 0.174
## SOAP_NDVI 172.771 291.392 0.593 0.568
##
## Residual standard error: 39.66 on 9 degrees of freedom
## Multiple R-squared: 0.342, Adjusted R-squared: 0.1958
## F-statistic: 2.339 on 2 and 9 DF, p-value: 0.1521
## Full model in respect to base canopy diameter
mod.full.bascalcanopydiam<-summary(lm(max.basalcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand))
(mod.full.bascalcanopydiam)
##
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI,
## data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.846 -23.404 -4.255 6.173 84.029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 856.50684 282.18923 3.035 0.0141 *
## SOAP_lidarCHM -0.08263 1.19195 -0.069 0.9462
## SOAP_NDVI -879.05838 340.38446 -2.583 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 46.32 on 9 degrees of freedom
## Multiple R-squared: 0.5097, Adjusted R-squared: 0.4008
## F-statistic: 4.678 on 2 and 9 DF, p-value: 0.04046
## Full model in respect to max canopy
mod.full.maxcanopydiam<-summary(lm(max.maxcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand))
(mod.full.maxcanopydiam)
##
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.166 -9.629 -3.622 0.498 49.964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.1373 116.7139 -0.130 0.900
## SOAP_lidarCHM 0.8595 0.4930 1.744 0.115
## SOAP_NDVI 4.3227 140.7836 0.031 0.976
##
## Residual standard error: 19.16 on 9 degrees of freedom
## Multiple R-squared: 0.3194, Adjusted R-squared: 0.1681
## F-statistic: 2.112 on 2 and 9 DF, p-value: 0.177
## Full model in respect to stem height
mod.full.stemheight<-summary(lm(max.stemheight~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand))
(mod.full.stemheight)
##
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.418 -15.603 -4.626 8.692 59.653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.6458 191.8097 -0.050 0.9610
## SOAP_lidarCHM 2.4095 0.8102 2.974 0.0156 *
## SOAP_NDVI -34.5459 231.3662 -0.149 0.8846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.49 on 9 degrees of freedom
## Multiple R-squared: 0.5603, Adjusted R-squared: 0.4626
## F-statistic: 5.735 on 2 and 9 DF, p-value: 0.02478
## Full model in respect to % score
mod.full.score.perc<-summary(lm(tree.score.perc~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand))
(mod.full.score.perc)
##
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22356 -0.06836 0.01661 0.09084 0.27052
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.446148 0.990811 1.460 0.1784
## SOAP_lidarCHM 0.008363 0.004185 1.998 0.0768 .
## SOAP_NDVI -1.248143 1.195144 -1.044 0.3236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1627 on 9 degrees of freedom
## Multiple R-squared: 0.3073, Adjusted R-squared: 0.1534
## F-statistic: 1.997 on 2 and 9 DF, p-value: 0.1916
